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ABSTRACT 

A two-dimensional flow governed by the incompressible Navier-Stokes equations with a 
steady spatially periodic forcing (known as the Kolmogorov flow) is numerically simulated. 
The behavior of the flow and its transition states as the Reynolds number Re varies is inves- 
tigated in detail, as well as a number of the flow features. A sequence of bifurcations is shown 
to take place in the flow as Re varied. Two main regimes of the flow have been observe : 
small and large scale structure regimes corresponding to different ranges of Re. Each of 
the regimes includes a number of quasiperiodic, chaotic and relaminarization windows. In 
addition, each range contains a chaotic window with non-ergodic chaotic attractors. Spa- 
tially disordered, but temporally steady states have been discovered in large scale structure 
regime. Features of the diverse cases are displayed in terms of the temporal power spectrum, 
Poincare sections and, where possible, Lyapunov exponents and Kaplan- Yorke dimension. 


'■Research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS1-18605 while the authors were in residence at the Institute for Computer Applications in 
and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1 Introduction 


It is now well established that a number of dissipative chaotic fluid flows 
may require relatively few dimensions for their description 1-4 . This has been 
examined perhaps in greatest detail for fluid models such as the Ginzburg- 
Landau (G-L) equation 5-8 , an equation which results from the study of the 
critical point in a variety of stability problems 9-12 . For closed fluid systems 
such as Rayleigh-Benard (R-B) convection and Taylor-Couette flow there is 
abundant evidence that the early stages of chaos are low dimensional. For 
R-B convection there also exists a numerical simulation which details the 
variation of dimension with Rayleigh number 13 . 

In this paper we consider a simple fluid flow which exhibits chaotic behav- 
ior and which is amenable to detailed examination. This flow, governed by 
the incompressible Navier-Stokes equations, is generated and maintained by 
a spatially periodic time independent forcing of the fluid. It is sometimes re- 
ferred to as the Kolmogorov flow (Arnold and Meshalkin 14 ). The stability of 
the primary flow has been reported on by Meshalkin and Sinai 15 and Green 16 . 
That flow is known to become unstable at a critical Reynolds number, Re c , of 
\/2 15,16 , beyond which a stationary cellular pattern appears 1 ' . At sufficiently 
high Re the cellular pattern itself becomes unstable and the resulting flow 
is unsteady and chaotic. Kraichnan 18 and Batchelor 19 have shown that two 
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dimensional turbulent flows exhibit an energy cascade toward low wavenum- 
bers (the infamous reverse cascade) and an enstrophy cascade toward high 
wavenumbers. This phenomenon is apparent in our calculations* An alter- 
nate view of this effect, in the present instance, is that the flow is unstable 
to long wavelength disturbances. The general instability of short wavelength 
structures to long wavelength disturbances has been investigated by Yakhot 
and Sivaskinsky, and under sufficiently anisotropic conditions, their study 
lead to the notion of negative viscosity for the long wavelengths 20,21 , since 
the growth rate is proportional to the square of the wavenumber. 

Our investigation is in large part numerical in nature. We examine in 
some detail the behavior of the flow and its transition states for a range of the 
bifurcation parameter Re. We use Poincare sections to examine the chaotic 
attractor of the flow, and apply the Kaplan- Yorke formula 22 to calculate 
the dimension of the chaotic attractor when it is possible. A number of 
remarkable features of the flow were discovered. These include windows of 
rclammarization in the parameter space, chaotic windows with Tion-cvgodic 
attractors and spatially disordered but temporally steady states. 

We will be caieful throughout this paper to use the term chaotic instead of 
the term turbulent . In particular we are referring to flows which when looked 
at in the time domain have broad band energy spectra. These flows also 
exhibit sensitive dependence on initial data (they have positive Lyapunov 



exponents). Thus by generally accepted criteria the flows are temporally 
chaotic. However, as will be seen, relatively few spatial Fourier modes are 
active and we say we have chaos rather than turbulence. Nevertheless the 
flow is chaotic in its Eulerian description and should not be confused with 
Lagrangian turbulence 23-25 in which the velocity fields are to a high degree 
laminar. 

Recently we became aware of the investigation of the Kolmogorov Flow 
by She 26 ’ 27 . While there is some unavoidable overlap, the two investigations 
are largely complimentary. She’s work is also mainly numerical and concen- 
trates on the high end of the Re range. (A minor difference is that he uses 
a higher value for the frequency of the forcing function.) The mathemati- 
cal tools used also differ. We will make extensive use of the power spectra, 
Poincare sections, Lyapunov spectra and flow visualizations for the investi- 
gation, while She relies mainly on time series and phase portraits. Several 
issues are viewed differently by us and we comment on these in the text. 
She’s further work done in collaboration with Nicolaenko is concentrated on 
the symmetry breaking homoclinic chaos 28,29 occurring beyond the Re range 
covered in this study. 
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2 Kolmogorov Flow 


We consider the forced incompressible Navier-Stokes equations in two di- 
mensions, 

ay + fy' = V' - u' = 0 (2.1) 

dJi u ' + Vp '/ P = v V'V + X sin £y'e x (2.2) 

In preparation for the normalization given below, dimensional quantities are 
denoted by primes, e x denotes the unit vector in the x-direction. This system 
has the solution or fixed point , 

u ' = U ' = sin ky'. (2.3) 

A Reynolds number for the flow may be naturally based on the maximum 
speed of this field, %/ vk^ and the length scale k ~ thus 

Re = */(fcV). (2.4) 

(A Reynolds could be based on the wavelength and maximum speed in which 
case 27 t Re is the Reynolds number). 

The lineal stability of (2.3) can be analyzed by writing 

u' = U' + 8u'(y)exp(iakx + at) (2.5) 

and investigating the resulting Orr-Sommerfeld equation gotten by lineariza- 
tion in the small perturbation, <5u' 30 . It is found that instability first sets in 
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( 2 . 6 ) 


at 16 ’ 16 

Re = Re c = \/2 

and is due to disturbances of infinitely long wavelength. Moreover, it can be 
shown that disturbances with ak > 1 never destabilize U'. 

Thus according to the stability analysis, the domain must be large in 
some sense in order for us to include growing modes. For computational 
purposes we fix the domain size and regard the spatial forcing frequency 
as the bifurcation parameter. With this in mind we choose the following 
normalization, 

u = i/k 2 u'/x , x = fcx'/n, t = xt'/{ n kv),p = pv 2 k 4 /x 2 , (2-7) 

where n is an integer which specifies the spatial frequency of the forcing. 
Under this normalization the Navier-Stokes equations become, 


du 

dt 


V * u = o 

1 


n 


+ u + — e x sin ny 




a 


with 


fl = nx/{vk) — nRe 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 


the bifurcation parameter. In this format the critical value of ft is 


ft r = nRe c = n\J 2. 


( 2 . 11 ) 
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Let us now briefly discuss the symmetry groups for the solutions of these 
equations which we will refer to as the Kolmogorov flows 31 . There are in fact 
two groups, one discrete and the other continuous. 

Referring back to the notation of (2.9), n refers to the number of vertical 
cycles in the horizontal forcing. If [u(x, y, <), v(x, y, t)] = u is a solution then 

g k u = [(-)*u ((-)**, y + k~^j , v ((-)**, y + *:£)] k = 0, . . . ,2n - 1, 

(2.12) 

are also (independent) solutions of the problem, i.e. they satisfy the equations 
and conditions of the problem. The generator of this cyclic group, g , is a 

glide-reflection 3 ~ of half wavelength. In addition rotation through tt is another 
group generator, i.e. 

Ru = [-u(-x, -y), -v(-x, -y)] (2.13) 

is also a solution. In all there are 4 n elements in the discrete symmetry 
group foi the Kolmogorov flow, thus for the case being considered, n = 4, 
we obtain a sixteen-fold increase in the number of possible solutions which 
we can view. 

The remaining symmetry group is the group of translations in the x- 
directions, i.e. if u is a solution then 

7',u = [u(.r + l,y),v(x + l,y)] (2.14) 

are also solutions to the problem for all 0 < / < 2tt. 
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3 Computational Method 


Based on the normalization given in the previous chapter we consider the 
computational domain 

0 < < 2tt (3*1) 

with periodic boundary conditions in the two directions (so the flow lies on 
the torus). 

The natural method of choice for the investigation of (2.8,9) under the 
27r-periodic boundary conditions is the pseudo-spectral method using Fourier 
expansions 33 . All derivatives are evaluated in the spectral space and updates 
are performed on the expansion coefficients of u(x,f). On the other hand the 
non-linear terms are evaluated in the physical space and then transformed 
back to the spectral space. The needed transformations are all efficiently 
done by the means of the Fast Fourier Transforms. 

In a variation on the usual formulation we have used the vector expansion 

set 

B„.(x)=(=^exp(.m.x), (3.2) 

I m I 

where m is the integer wavenumber pair (??i 1 ,m 2 ). Using such a basis we can 
expand the velocity field u as 

u(x,t) = ^u lu (/)B m (x) (3.3) 

111 
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where the sum is performed over all the components of the wave vector m. 
The basis elements are easily seen to be divergence free, 


v * B m = 0 (3.4) 

and complete in the space of square integrable 27r-periodic divergence-free 
vector functions. An advantage of this formulation is that gradients of scalar 
functions (such as the pressure term in the Navier-Stokes equations) are 
orthogonal to the divergence-free space 34 and so a result of a Galerkin for- 
mulation is that pressure is eliminated entirely from the problem. A compu- 
tational method of this type was used be Leonard and Wray 35 in their study 
of the stability of pipe flows. 

Throughout our study the underlying forcing frequency was taken to be 
n 4 and so that the critical value of the bifurcation parameter Q, is 


Q c = 4v/2. 


(3.5) 


Furthermore, we shall be interested in the long-term, asymptotic behavior of 
the flows developing from the primary flow 


TT sin iy 

u = u = -?r e ~ 


(3.6) 


which is a fixed point solution of (2.8,9). To facilitate this study we use a 
small fixed perturbation of the primary flow including all positive and nega- 
tive wavenumbers as the initial conditions for the calculations. In particular, 

S 


we set 


u 0 (x, 0) = ^ sin 4ye x + ^ [(l + 0 B m(x) (1 i)B_ m (x)] (3.7) 


mi >0 
m 2 >0 


If we write 


u 0 = ( u 0 ,v 0 ) 


(3.8) 


then for example u 0 has the form 

l 


1 

Uo = n 


V 


mi >0 
m 2 ^0 


m 2 

m 


sin 4y + 10 -4 t— M caKmi® + m 2 y) - sin^x + m 2 y) 


(3.9) 


The computations were performed for the range 


a/n c < i2.5 


(3.10) 


After allowing a considerable time for the calculation to settle down (on 
average 5000 time units), we examine the resulting flow on the attractor. It 
was found out that a 32 x 32 spatial grid was adequate for the resolution 
at the ratios (3.10). Two types of dealiasing were used during the course of 
the simulation as were grids of twice the resolution. Neither were found to 
significantly alter the results. 
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4 Tools of the Analysis 

Here we briefly comment on the tools used in the analysis of the di- 
rect numerical simulation of the Kolmogorov flow. The tools needed for the 
investigation can be grouped according to their function. 

4.1 Temporal Power Spectra 

We need some mathematical tools to distinguish between chaotic and 
quasiperiodic regimes of the flow. The temporal power spectra (PS) are 
indispensable for this purpose 36 . 

Let PS(j(t)) denote the temporal PS of the time signal 7 (f). It follows 
from equation (3.3) that there are 32 2 = 1024 individual time histories avail- 
able for the use in the temporal PS. For the problem at hand, we define 
the energy PS box averaged per mj wavenumber in x-direction as (modulo 
a multiplicative normalization) 

PS( u) (mi) = ^2 /' > 5(/?e(u( mi m2 )(<))) -f PS(Im(u {miiTn2) (t))). (4.1) 

m2 

In other words, we take the temporal PS of all the time histories result- 
ing from the projections onto basis function ( 3 . 2 ) which include irq in the 
wavenumber pair m = (mi,m 2 ) and average over them. This is in keeping 
with the definition of the PS of a time stationary signal as being the Fourier 
transform of the autocorrelation. In the calculations we frequently set nq to 
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zero and use the resulting temporal power spectrum to distinguish between 

chaotic and quasiperiodic regimes of the flow. 

For later reference we remark now that since the PS is the Fourier trans- 
form of the autocorrelation it follows that as the frequency | / |— » oo the PS 
vanishes faster than any inverse power of | / |. This is a simple consequence 
of the fact that in any numerical simulation we are integrating a system of 
ordinary differential equations with smooth direction fields and hence the 
solution is C°° in the time variable. 

4.2 Poincare Sections 

We use Poincare sections to study geometrical structure of the attractor 
and to qualitatively discuss its degree of complication. Methods for the 
calculation of these sections have been given by Keefe 5 and others. The 
phase space employed in the calculation of the Poincare sections is spanned 
by the spatial Fourier modes of the solution u(x,f) 

u(x,i) = y^Um(QB m (x) (4-2) 

III 

where u m and B m are complex. 

The choice of the sectioning plane for the Kolmogorov flow becomes very 
important, for it determines which attractor features are covered by the 
study. Here, we arc mainly concern with the visual structure of the at- 
tractor itself. CPU expenditures to accumulate', a large enough collection of 
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points for a highly resolved Poincare section become considerable. Therefore 
a need arises for the use of the discrete symmetry group of the flow to extend 
the available data. This in turn suggests that some care should be exercised 
in the choice of a sectioning plane. In particular, an acceptable choice is a 
plane which is invariant under the symmetry group. It is also a good practice 
to choose as a section, a plane for which the flow projections maximize the 
energy. With that in mind, we take the forcing function plane 

/m(u ( o, 4) ) = D (4.3) 

/m(u (0t4 )) > 0 ( 4 . 4 ) 

as our invariant under the action of the symmetry group section plane for the 
Poincare sections used later. The choice of the constant D controls apparent 
positioning of the hyperplane in the phase space and is dependent upon Q. 
In actual calculations, we consider a series [/m(u (0t4) (<,-))]£, of snapshots of 
the variable Im{u (0A] (t)). We then set D to its average, i.e. 

1 N 

D = Jr 2ZM«(o,4) (*.))• (4.5) 

» = 1 

The Poincare sections shown in the later sections plot Re(u( 0 jj) vs /m(u (0 ,*). 

For further reference, we note here representation of the symmetry group 
(2.12,13) in terms of the spectral components of u(x, t). If u(x, t) is rep- 
resented as m (4.2) then the glide-reflection $tu(x,*) (2.12) can be written 
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< 7 ll(x,t) = U (-n»j,m 2 )(0 e ' " 7r B m (x) = 5Z U (”M.- m 2)^^ e " ® m ( X ) 

v ^ ' m 

111 111 

(4.6) 


while the rotation Ru (2.13) is 


i?u(x, t) = ^2 u -m(0Bm( x ) — u m(0®m( x ) (4-7) 

m m 

It then follows from (4.5) and (4.6-7) that our sectioning plane is invariant 

under the symmetry group of the flow. 

From the examination of the equation (4.7), it is apparent that the ro- 
tation corresponds to the reflection around imaginary axis for any variable 
Um (t). On the other hand glide-reflection action (4.6) of the symmetry group 

on tt( O ,i)(0 gives 


- «(o,i )(0 e ' 4 

This corresponds to fx rotation in i2e(«( O ,i)(0) ~ /™(«(o,i)(0) P lane - rlius 
the discrete symmetry group of the flow induces an extra symmetries onto 
U ( 01 )(t) plane. They correspond to the reflection around 1 g x, i = 0, 1,2,3 
axis. These symmetries figure prominently in the discussion of the results of 
our numerical simulation. 
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4.3 Lyapunov Spectrum 


All of the mathematical tools mentioned above can indicate the presence 
of chaos in the dynamical system, but do not directly address the level of 
randomness present, or how chaotic a solution is. 

We compute the Lyapunov exponents for the two dimensional Navier- 
Stokes equations (2.8,9) to study the dimension of the attracting set of the 
chaotic flow. We write the governing system in symbolic form as 

du -r-/ \ 

Tt = r(u) - < 4 - 9 ) 

and denote a reference or fiducial trajectory by u°. If an initially nearby 
trajectories is written as 


u = u° + ev 


(4.10) 


then in the limit e — » 0, v satisfies the variational equation, 

d 


-v = £(u> 


where 


£(u“) 


<5u° 


(4.11) 


(4.12) 


is easily obtained through the linearization. 

1 he laigcst, Lyapunov exponent is defined as the maximum exponential 
i ate at which nearby trajectories diverge 


A 



(4.13) 
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Thus the idea of a Lyapunov exponent can be regarded as a generalization 
of the stability exponent of an equilibrium solution, and finding a positive 
Lyapunov exponent is an unambiguous signature of a chaos. 

In a practical sense we expect all initial data, v(0), to lead to A 1? i.e., 
as time evolves v will become parallel to the direction of greatest growth. 
Thus to obtain the next greatest growth rate, and so forth can be a touchy 
problem. Algorithms for calculating these has been given by Bennetin et 
al 37 , Shimeda and Nagashima 38 and Wolf et al 39 , and are further explained 
by Fitzmaurice 40 . The essence of the trick lies in computing volumes of 
infinitesimal ellipsoids. In fact if c n (r) denotes the volume of an n-ellipsoid 
it then follows that 

A x T A 2 T . . . T A n = lim In c„(t). (4-14) 

T — ►OO 7- 

Thus after (4.13) one calculates Ax + A 2 , then Ax + A 2 + A 3 , etc and thereby 
the spectrum Ax,A 2 , — 

Since the system in question is dissipative we know that the sum (4.14) 
will become negative for n large enough. And since it is chaotic it must be 
positive for n = 1. Therefore this sum must cross zero at some point as one 
increases u. One may then use the formula of Kaplan and \ 01 ke 

d,, = N - (4-15) 

A/v+i 
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to estimate the dimension, d £, of the attractor where N is the largest integer 
for which the sum is positive. Constantin et al 41 have shown that the Kaplan- 
Yorke formula gives an upper bound for the Hausdorf dimension. 
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5 Results 


In this section we present the results of our direct numerical simulation 
of the Kolmogorov flow. We have found that there are two main regimes 
of the flow corresponding to different ranges of the bifurcation parameter 
fl. Tentatively these are termed as being the regimes of small and large 
scale structures. Explanation for this terminology will be evident from the 
streamline plots to be discussed. 

5.1 Small Scale Structures 

0 < SI/ Sic < 1 

For this range the asymptotic solution of the N-S equations (2.8,9) is 
the fixed point (3.6). Streamlines and vorticity plots for the parallel shear 
flow in horizontal direction are shown in Figure 1. 

1 < SI/ Sl c < 1.97 

As we increase 0, the horizontal shear flow becomes unstable and the 
flow pattern changes. At Sl/Sl c = 1, the parallel shear flow bifurcates to 
a new solution. In this range the motion is steady and purely cellular and 
illustrations of this are shown in Figures 2 and 3. It is interesting to note that 
at the lower range of the parameter f l the flow is almost horizontal (Figure 
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2a) while at higher values of 0 in this range the motion is mainly in the 
vertical direction (Figure 3a) with well defined cells. In the first instance the 
saddle points of the separatrices are formed horizontally while in the latter 
case they are joined vertically. This is in an agreement with the experiments 
of Bondarenko, et al 42 for a magnetohydrodynamic problem under similar 
forcing. 

1.97 < n/n c < 2.22 

In this range the flow still have the steady cellular structure but the 
number of cells in the horizontal direction have doubled. The new flow 
patterns are depicted in Figure 4. Explanation of this phenomena relies on 
the fact that the horizontal wavenumber of the fastest growing mode changes 
from 1 to 2, also reported by Green 16 . A very important feature of the flow 
in this and previous ranges is that it is invariant under the discrete symmetry 
group transformations. By invariance we mean that for any transformation 
G in the discrete symmetry group of the flow there is 0 < / < 2tt such that 

T,Gu = u (5.1) 

where 7/ is a translation by l in x-direction and u is a solution of N-S equa- 
tions (2.8,9). In words, we say that the flow is invariant, under the discrete 
symmetry group if, after application of any transformation from the discrete 
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symmetry group to the flow, we can translate the transformed flow in hon 
zontal direction in such a way that it coincides with the original flow. 

2.2 < 0/0 c < 3.53 

As $7 increases, small amplitude perturbations at the saddle points 
will lead to the destabilization and eventual breakdown of the steady cellular 
pattern. At Sl/Sl e = 2.2 another bifurcation in parameter space takes place. 
For the range under discussion there is a gradual breakdown in the symmetry 
of the flow. The new flow patterns are depicted in Figure 5. It is perhaps 
remarkable that these cellular structures are still steady. 

3.53 < ft/ftc < 3.6 

At Q/n c = 3.53 the statistically stationary flow is no longer steady. 
Now, the flow spends most of the time in the steady cellular state shown 
in Figure 5. But from time to time, the flow is intermittently disturbed by 
a random horizontal shifts occuring on a very rapid time scale. The time 
between bursts is very long and is a decreasing function as the parameter 
S7 reaches the upper end of the range. Figure 6 shows the time series of 
Re(u( o,i)). Here, the plateaus correspond to a steady cellular state. After the 
horizontal shift, the cells oscillate with a small amplitude foi a short while, 
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and then settle down very rapidly to a steady pattern before undergoing the 
next horizontal shift. This state is chaotic but only over a very long time 
scale, i.e. the Lyapunov spectrum has some very small positive exponent. 

3.6 < n/n c < 4.47 

In this range the nature of the flow changes once again. We now find a 
periodic oscillations occuring in the cellular structures as depicted in Figure 
5. Also, the solution develops a slow drift in horizontal direction, but without 
a net transport off mass. The drift may be regarded as the wave motion. It 
also follows from symmetry that the drift wave can move in either direction. 
We conclude from f igure 7 which shows the temporal power spectrum that 
there is a principal frequency, the large peak, and subharmonic frequencies 
of it present in the flow. These correspond to the oscillations of the pattern 
around larger vortices and inner oscillations of the vortices themselves. It is 
interesting to note that we have found no regime of the flow when only the 
principal frequency is present. 

4.47 < fl/n c < 4.58 

Now, the flow enters its second chaotic regime. Figure 8 depicts a 
typical time series of the signal Re(u {0 ^) in this range, from which we con- 
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elude that the chaos is intermittent. Examining the power spectrum shown 
in Figure 9, we see that the principal and subharmonic frequencies are promi- 
nently present in what is a broad band spectrum. We did not attempt to 
calculate the Kaplan- Yorke dimension for the dynamical system in this range 
of ft because of the very long time interval between the bursts of activity in 
the flow. A most remarkable feature of the flow in this range is that the 
chaotic attractor is not ergodicl The dynamical system has several chaotic 
attractors each with its own basin of attraction lying quite close to the fixed 
point (3.6). On the basis of extremely long calculations we conclude that 
these attractors are not connected. Figure 10 shows the Poincare section 
depicting cuts of the chaotic attractors. Two darker points in the blown up 
section Figure 10b correspond to the punctures produced by the solution in 
the laminar phase of the intermittency. In actual calculations we obtain only 
one of the attractors (depending on the chosen initial data) , and the rest 
are obtained by using the discrete symmetry group. Alternatively, we can 
apply the symmetry group to the initial data and thus generate each of the 
attractors. Another peculiarity of the flow is that there are only eight sym- 
metric attractors instead of sixteen as predicted by the discrete symmetry 
group. This happens because the action of five glide- reflections followed by 
the rotation by 7 r applied to the solution maps the attractor onto itself. This 
induces an extra symmetry into the attractor corresponding to the reflection 
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around the axis denoted by the dashed line in Figure 10b, as it is discussed 
in the earlier part of this paper dealing with the Poincare sections. 

4. 58 < n/n c < 4.65 

As Q increase, different symmetric attractors start intersecting each 
other and now the flow is intermittenly chaotic with one ergodic attractor 
consisting of eight connected symmetric components. Figure 11 shows the 
time series of i?e(u( 01 )). The evidence of the intermittent bursts is clearly 
present and sudden jumps correspond to the solution jumping to a different 
component of the attractor. We present the temporal power spectrum in 
Figure 12. This new state is more chaotic than the previous case, but still 
there is evidence of the principal frequency from the lower parameter values. 
The Poincare section is shown in Figure 13 (“butterflies” section). Darker 
regions in the blown up section of the attractor in Figure 13b correspond to 
the points where the solution spends most of the time between intermittent 
bursts. Moreover, comparing Figures 13 and 10, we notice that the main 
features of the non-ergodic attractor are still present. 

4.65 < n/n c < 4.71 


Another remarkable feature of the Kolmogorov flow is, that counter 



to intuition, the parameter space fi contains windows of relaminarization 
between different chaotic ranges of the flow. This range is an example of 
one of such windows. Here, the flow returns to a periodic regime with a 
principal and subharmonic frequencies in the power spectrum. Again, there 
is no regime with just one frequency. This flow is similar to the range 3.6 < 
Q,/Q, c < 4.47. It also, develops a slow steady drift in the horizontal direction. 
We also note that this drift appears to vanish in the chaotic windows of the 
flow. 

At ft = 4.72ft c we see the formation of large scale structures. 

5.2 Large Scale Structures 

Although we are mainly interested in the asymptotic behavior of the 
Kolmogorov flow, it is of interest to examine how the large scale structures 
are formed. The process is basically the same for all values ft > 4.72ft c . 

We start with a small perturbation of the fixed point flow (3.6). This state 
persists for roughly the time scale predicted by the linear stability analysis. 
Then on a very rapid time scale the flow changes to the unsteady cellular 
pattern of the sort illustrated in Figure 5. This state with many small cells 
persists for a very short time and then the flow plunges into chaos. 1 he new 
state is vigorously chaotic, characterized by the rapid movements of a small 
vortical structures, a snapshot of which is shown in Figure 14. By a rough 



count it is seen that the same number of vortical structures appears in Figure 
5 and 14. We refer to this state as metastable or transient chaos following the 
terminology of Grebogi, Ott and Yorke 43 . Of all the chaotic states observed 
in our numerical experiments, this is the only state which can be called 
chaotic both spatially and temporally and hence is turbulent according to 
common convention. At the lower range of ft this metastable chaos holds for 
a relatively long time, while at a higher values of the bifurcation parameter 
ft it can be extremely short. 

In any case the phase point finds an escape hole in the metastable attrac- 
tor and we eventually see the formation of a large scale structures such as 
those depicted in the Figure 15. The new attractor is stable. Considerable 
vortex merging has taken place in passing from Figure 14 to Figure 15. The 
two large structures shown in the streamline plot of Figure 15a are some- 
what misleading. On a closer examination of the vorticity plot, we see that 
we have roughly the same number of distinct vortices as before and that the 
large structures in the streamline plot consist of a number of a counterrotat- 
ing vortices. These new large scale structures are completely different from 
the small scale structures which we dealt in the previous subsection. Exam- 
ining the vorticity plot shown in Figure 5b, we notice that at the small scales 
the vortices arranged in a lattice. Thus, the unsteady small scale structures 
flow can be considered as oscillations or “breathing” of the lattice structure. 
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In the metastable regime this lattice structure is broken, and we have a soup 
of vortices moving randomly around. Later the vortices merge together and 
arrange themselves on a curved line with two vortices containing most of the 
circulation. These vortices correspond to the eyes of the large scale struc- 
tures depicted in the streamline plot of Figure 15a. It is of interest to note 
here that the vertical position of the eyes of the large scale structures are at 
hyperbolic points of the forcing functions when the flow is not chaotic. 

We will now describe in detail the bifurcations occuring in the large scale 
structure regime of the flow. 

4.72 < ft/ft c < 5.9 

In this range the flow is periodic with one frequency. In addition, there 
is a horizontal drift present in the solution as also noted by She 26 . As ft/ft c 
approaches 5.9 both the frequency of oscillations and drift speed approach 
zero. Figure 16 shows the temporal power spectrum. 

ft = 5.9ft c 

At ft = 5.9ft c , the flow reaches the steady, spatially disorganized state 
depicted in Figure 15. The flow is chaotic only in the spatial domain! 
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5.9 < n/n c < 9.47 


As we start increasing ft further, the flow develops both periodic oscil- 
lations and a slow drift in the horizontal direction. The direction of the drift 
changes according to the particular attractor the flow settles in. Moreover, 
as Q increases, so does the drift speed and frequency of the oscillations. The 
periodic oscillations can be described as the “breathing” in the inner shells 
of the large scale structures. At S7 = 9.042fl c a period doubling bifurcation 
takes place in the oscillations of the inner structures. Now the principal and 
subharmonic frequencies correspond to the breathing of the large scale struc- 
ture itself and the oscillations of the inner shells. This regime persists up 
to = 9.45ft c , where another period doubling bifurcation in the oscillations 
takes place. We can see clearly the principal and subharmonic frequences in 
Figure 17 showing the temporal power spectrum of the flow. 

9.47 < fl/Q c < 9.66 

At fl = 9.47H C the flow enters another chaotic regime which continues 
until Q reaches 9.66fi c . We again find the remarkable feature in this win- 
dow that there is chaos but that the attractors are non-ergodicl Here, the 
chaos is present in the irregular breathing-like oscillations of the inner struc- 
tures which induces random oscillations of large scale structures. Wc observe 
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in Figure 18, which shows the temporal power spectrum, that the driving 
frequency and its subharmonics still dominate the motion. Horizontal drift 
diminishes to a trickle and sometimes changes directions. On average, there 
is no drift at all. Figure 19 depicts the Poincare section of the attractor 
showing eight symmetric disconnected attractors. The slopes of the lines in 
Figure 20 correspond to the values of Lyapunov exponents and the calculated 
Kaplan- Yorke dimension is 7.9. 

9.66 < f l/Cl c < 10.805 

At Cl = 9.66H C the flow enters another laminar window in the pa- 
rameter space Cl. Here, the solution slowly drifts in the horizontal direction 
and has the principal and subharmonic frequencies. At Cl/Cl c = 10.4 a pe- 
riod doubling bifurcation takes place. Figure 21 shows the temporal power 
spectrum at Cl — 10.7Q C . 

10.805 < Cl/Cl c < 12.5 

As we increase Cl even further, the flow enters another chaotic regime at 
Cl = 10.S05f! c . The flow now appears to be governed by an ergodic attractor 
which consists of eight symmetrical components and the solution intermit- 
tently jumps from one component to another. There is a very long time 


delay between successive jumps at ft slightly greater than the 10.805, but 
it rapidly decreases as ft increases. The driving frequency and its subhar- 
monics dominate the power spectrum shown in Figure 22. An examination 
of Figure 23 showing the Poincare section of the chaotic attractor, indicates 
that the components of the attractor have a fractal structure. The calculated 
Kaplan- Yorke dimension is 8.9. Another interesting point is that the direc- 
tion of the drift changes with each successive jump to a different component 
of the attractor. 

As we continue to increase ft, the different components of the chaotic 
attractor merge together and become indistinguishable. Nevertheless, the 
driving and subharmonic frequencies are still dominate the temporal power 
spectrum shown in Figure 24 for ft = 12.5ft c . Now, the Kaplan- Yorke di- 
mension increases to 10.8. Thus we might hope to describe the system in 
terms of a dynamical system having many fewer than roughly 1000 degrees 
of freedom that we have been allowing in our spectral simulations. 

As in the case of the small scale structures, the average drift speed di- 
minishes to zero for the chaotic regimes and the drift is not present at all for 
ft > 12.5ft c . 

This concludes oui detailed description of the behavior of the fixed point 
(3.6) of the Kolmogorov flow as a function of the bifurcation parameter ft. 
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6 Summary 


In this paper we investigated in some detail a simple fluid flow governed by 
the incompressible two-dimensional Navier-Stokes equations with a spatially 
periodic time independent forcing. Two main regimes of the flow have been 
observed: small and large scale structure regimes corresponding to different 
ranges of Reynolds number Re. A sequence of bifurcations takes place in 
each of the regimes of the flow. 

The first three bifurcations taking place in the small scale structure regime 
are steady: the first corresponds to the appearance of the cellular structures, 
the second doubles the horizontal wavenumber of the cellular structures, 
while the third breaks the symmetry of the flow. In addition, three chaotic 
windows are discovered in the small scale structures regime. The first window 
is chaotic over a very long time scale. The second window is an intermittency 
with a remarkable feature that the chaotic flow has non-ergodic attractors. 
Finally, the last chaotic window shows intermittency with one connected 
chaotic attractor consisting of eight symmetric components (“butterflies”). 
There are laminar windows between the first and the second chaotic regimes, 
and after the last chaotic regime. 

The large scale structure regime contains a number of quasi-periodic and 
chaotic windows. Before the formation of the large scale structures, the flow 
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spends considerable time in metastable chaos which approaches true turbulent 
flow for the range of Re under the discussion. A spatially disordered, but 
temporally steady state was observed between first two laminar windows. 
Two period-doublings take place in the second laminar window. The chaotic 
regime which follows has the same remarkable property of chaotic but non- 
ergodic attractors as is the case with the small scale structure regime. In 
addition, there is a window of relaminarization with a quasi-periodic flow 
between the two chaotic regimes. 

In the future, we plan to use the proper-orthogonal decomposition 44 to 
extract the coherent structures of the flow and study the resulting low order 
dynamical systems 45 . 
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1. Parallel shear layer flow corresponding to the fixed point solution at 
fi = 0.5Q C . a) streamlines. 
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1. Parallel shear layer flow corresponding to the fixed point solution at 
fl = 0.5fi c . b) vorticity. 
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2. Steady cellular flow found at Q = 1.12ft c with the motion mainly in 
horizontal direction, a) streamlines. 
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2. Steady cellular flow found at fl = 1.120= with the motion mainly in 
horizontal direction, b) vorticity. 











3. Steady cellular flow found at Q = 1.5fi c with the motion mainly in ver- 
tical direction and the fastest growing perturbation with wavenumber 
1 in horizontal direction, a) streamlines. 


42 


ORIGINAL PAGE IS 

OF POOR QUALITY 




Steady cellular flow found at Cl = 1.5f2 c with the motion mainly in ver- 
tical direction and the fastest growing perturbation with wavenumber 
1 in horizontal direction. b)vorticity. 
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4. Steady cellular flow found at Q = 2fl c with the motion mainly in ver- 
tical direction and the fastest growing perturbation with wavenumber 
2 in horizontal direction, a) streamlines. 
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4. Steady cellular flow found at ft = 2Q C with the motion mainly in ver- 
tical direction and the fastest growing perturbation with wavenumber 
2 in horizontal direction, b) vorticity. 
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5. Steady cellular flow found at ft = 3Q C . a) streamlines. 
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5. Steady cellular flow found at ft = 3J2 C . b) vorticity. 
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6. Time series of i2e(u( 0| i))(i) at fi = 3.56Q C . Plateaus correspond to a 
steady cellular flow. 
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7. Power spectrum of the velocity u at ft = 4ft e . The principal frequency 
and its subharmonic dominates the spectrum. 
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8. b) blow up of a). 


51 





Log Power 



Frequency / Hz 


9. Power spectrum of the velocity u at (1 = 4.48Q e . The flow is chaotic, 
but the principal frequency and its subharmonics still dominate the 
spectrum. 
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11. a) Time series of fle(u (0 ,i ))(i) at fi = 4.63ft c . The flow is intermittent 
chaotic. The sudden jumps correspond to the solution jumping from 
one component of the chaotic attractor to another. 
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13. b) One of the components of the attractor. 


59 




14. Ins t an taneous flow pattern found in met els table chaotic regime at fl = 
10fi c . b) vorticity. The flow is both spatially and temporally chaotic, 
and thus may be called turbulent. 
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15. Steady spatially chaotic flow found at ft = 5.9fl c - b) vorticity. Also, 
this is an example of laxge scale structures. 
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18. Power spectrum of the velocity u at ft = 9.47fl c . The flow is chaotic, 
but the driving frequency and its subharmonics are still dominant in 
the spectrum. 
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19. a) Poincare sections of the chaotic attractor at Cl = 9.47fl c . Again, 
attractors are not ergodic. 
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22. Power spectrum of the velocity u at fi = 10.80512 C . The flow is chaotic, 
but the driving frequency and its subharmonics are dominant in the 
spectrum. 
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23. a) Poincare section of the chaotic attractor at fi = 10.80517, 
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23. b) Blow up of one of the components of the attractor. The fractal 
structure of the attractor is evident from blow up. 
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24. Power spectrum of the velocity u at fi = 12.5ft c - The driving frequency 
and its subharmonics Me still dominant in the spectrum. 
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